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Abstract 

Background: A major challenge in mathematical modeling of biological systems is to determine how model 
parameters contribute to systems dynamics. As biological processes are often complex in nature, it is desirable to 
address this issue using a systematic approach. Here, we propose a simple methodology that first performs an 
enrichment test to find patterns in the values of globally profiled kinetic parameters with which a model can 
produce the required system dynamics; this is then followed by a statistical test to elucidate the association 
between individual parameters and different parts of the system's dynamics. 

Results: We demonstrate our methodology on a prototype biological system of perfect adaptation dynamics, 
namely the chemotaxis model for Escherichia coli. Our results agreed well with those derived from experimental 
data and theoretical studies in the literature. Using this model system, we showed that there are motifs in kinetic 
parameters and that these motifs are governed by constraints of the specified system dynamics. 

Conclusions: A systematic approach based on enrichment statistical tests has been developed to elucidate the 
relationships between model parameters and the roles they play in affecting system dynamics of a prototype 
biological network. The proposed approach is generally applicable and therefore can find wide use in systems 
biology modeling research. 
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Background 

Systems biology aims to unravel the design principles 
of living organisms from a network perspective [1,2]. 
Advances in experimental studies have generated a 
large amount of data on several key biological processes 
[3-8], and networks of interactions between molecular 
species have been hypothesized [9-14]. Despite these 
advances, one unresolved challenge in systems biology is to 
understand how the hypothesized molecular interactions 
can lead to the observed biological phenomenon for com- 
plex biological systems. One way of pursuing this is via 
mathematical modeling of biological processes, which can 
also generate testable hypothesis for future experiments. 
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Biological processes are often complicated, and the 
complexity of their mathematical models usually in- 
creases with the amount of parameters involved. This 
generally gives rise to two fundamental problems in 
mathematical modeling. First, it is possible to have mul- 
tiple sets of parameter values that are equally likely to 
produce the observed data and finding the "best" param- 
eter set might be insufficient to fully characterize a bio- 
logical system if such a parameter set is not the only set 
with biological relevance [15]. Second, understanding the 
role of individual parameters on different aspects of ob- 
served systems dynamics can be difficult for parameter- 
rich models as it might be too time consuming to explore 
this systematically and exhaustively. Although a number of 
approaches, for instance the genetic algorithm-based method 
[16] and others [15,17], have been developed to search for 
parameter solutions in high dimensional spaces, they have 
not been extended to make inferences on the contribution of 
individual parameters to specific components of the system 
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dynamics. As such, in this paper we propose a general frame- 
work for addressing the aforementioned problems in math- 
ematical modeling of biological systems. Our methodology 
can be summarized as follows. For a given mathematical 
model of a biological process, one first defines the systems 
dynamics that the model is required to reproduce; one then 
searches through the parameter space by a sampling method, 
and keeps those sets of parameters with which the model 
can produce the required dynamics. These then form a 
matrix Z of functional parameter values, with each element 
Zy denoting the value of the ;-th parameter in the z-th param- 
eter set. This matrix then undergoes statistical analysis to test 
whether a particular parameter is biased towards a certain 
value (or certain range of values) for the model to produce 
the target dynamics. After this is done for all parameters, 
the results can be compiled to identify recurrent par- 
ameter values and any patterns they might form. Fi- 
nally, if the systems dynamics can be decomposed into 
different components or parts, then further analysis can 
be performed to associate a parameter with a particular 
aspect of the dynamics. 

Our methodology is demonstrated on a well-studied 
example of adaptation dynamics, that of the chemotaxis 
of Escherichia coli [18-21]. Generally speaking, adapta- 
tion refers to the ability of an organism adjusting to a 
new environment, and it is thought to be an important 
attribute for survival under fluctuating conditions 
[22-28]. In Escherichia coli, adaptation allows its 
chemotaxis system to reset stimulus-induced output to 
pre-stimulus value, even upon persistent external 
stimulation [29,30]. The dynamics of chemotaxis adap- 
tation has two parts (Figure 1A): first, the output signal 
of the system exhibits a sharp increase after the initial 
stimulation, and this is referred to as the sensitivity 
phase; second, after the initial sharp rise, the output 
signal decays to its initial state, and this is referred to as 
the precision phase. Figure IB illustrates the molecular 



processes involved, which have been identified experimen- 
tally [31,32]. Briefly, input signals are fed into the histidine 
kinase CheA-chemoreceptor R complex, and CheA then 
phosphorylates CheY (which can be dephosphorylated by 
CheZ) to regulate the process that drives the flagella of a 
bacteria. The key point here is that the activity of CheA is 
determined by the level of methylation of the CheA- 
bound receptor complex, which is controlled by demethy- 
lase CheB, which is in turn positively regulated by CheA 
itself, forming a negative feedback loop (Figure IB). Re- 
cently, Ma et al. [33] constructed a mathematical model 
for a three-node enzyme network and found that only two 
major types of network topologies can produce dynamics 
associated with adaptation (Figure 1C and Additional file 
1: Figure S8). One topology consists of a negative feed- 
back loop with a buffer node (NFBLB for short): node A 
positively influences the production of nodes B and C, and 
node B in return negatively regulates node A (Figure 1C). 
The other topology has an incoherent feed-forward loop 
with a proportioner node (IFFLP for short): here, node A 
induces node B, which in turn induces node C, and nodes 
A, B and C also have inhibiting role on nodes C, A and B, 
respectively (Additional file 1: Figure S8). The enzyme net- 
work driving the chemotaxis of E. coli has been found to 
resemble the NFBLB model [33]. Thus, in the rest of this 
paper, we use the NFBLB model to demonstrate our 
methodology in order to better understand how individual 
parameters contribute to the mechanism underlying the 
chemotaxis of E. coli Empirical findings from the litera- 
ture will be compared to our numerical results for valid- 
ation purposes. Results for the IFFLP model will be 
presented in supplemental data. 

Methods 

Model of the adaptive enzyme network 

The original E. coli chemotaxis model was proposed as a 
two-state model [29] and was later expanded to include 
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Figure 1 Dynamics and models of an adaptive enzyme network. (A) A schematic illustration of adaptation dynamics: sensitivity refers to the 
magnitude of the change in the output after the introduction of an external signal, and precision refers to the ability of the system returning to 
its pre-stimulus state after being perturbed by the external signal. Two quantities, a sensitivity score and a precision score, can be defined to 
measure these two dynamic properties (see Methods for details). (B) Network topology of the chemotaxis machinery in E. coli. (C) An enzyme 
network with a negative feed-back loop, known as the NFBLB model, that exhibits similar topology to the E. coli chemotaxis circuit where 
v n i W n2 ) (n = A, B or Q represents the activation (deactivation) process of the rate equation for node n. 
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the phosphorylation cascade [34]. In this study, the 
model we used is essentially the same as that used by [29] 
and [34], in which we consider the CheA-bound receptor 
complex as a single entity (node A in Figure 1C) and as- 
sume that this receptor complex exists in either a CheA 
active (A m ) or a CheA inactive (A) state. The superscript 
m denotes the methylated form of the receptor complex. 

The binding kinetic equation for the inactive receptor 
complex is given by, 
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where / denotes the concentration of chemo-attractant, 
Ki A is the ligand dissociation constant and k IA is the lig- 
and catalytic constant. 

The demethylation kinetic equation for the active re- 
ceptor complex is given by, 



B P + A mJ^ B P A m B P + a 



(2) 



where B denotes the concentration of phosphorylated 
demethylase CheB, K BA is the dissociation constant of 
phosphorylated CheB and k BA is the catalytic constant of 
phosphorylated CheB. 

The process of phosphate group transfer from active 
CheA to CheB is given by, 



(3) 



where K AB and J<ab are the dissociation constant and the 
catalytic constant of active CheA, respectively. 

The dephosphorylation of phosphorylated CheB is 
given by, 



B p ^B, 



(4) 



where k FB B is the dephosphorylation constant. 

Transfer of the phosphate group from active CheA to 
CheY (node C in Figure 1C) is given by, 



(7) 

where / denotes the input signal (i.e. the concentration 
of chemo-attractant); A, B and C denote the concentra- 
tions in active states; (l-A), (1-B) and (1-C) denote the 
concentrations in inactive forms; F B and F c are the con- 
centration of deactivating enzymes (assumed to be a 
constant value of 0.5 as in [33]) that transform the active 
states of B and C into inactive states. Kinetic parameters 
and are respectively the catalytic rate constant and 
Michaelis-Menten constant for the catalytic reaction be- 
tween substrate i and its regulator (activating or deacti- 
vating) enzyme /, where i and j = A, B, C, F B , or F c . Note 
that for each node A } B or C, the first term (vj) and the 
second term (v 2 ) of the equation represent its activation 
and deactivation rates respectively. 

Measuring adaptation 

Two quantities were used to evaluate the performance 
of a kinetic parameter set in producing adaptation dy- 
namics: (i) sensitivity to the input stimulus (equation 
(8)), which is defined as the difference between output 
response and the initial steady-state value, and (ii) preci- 
sion (equation (9)), which is defined as the inverse of the 
difference between pre- and post-stimulus steady state 
values. The corresponding mathematical equations for 
these two quantities are [33]: 



A m + c J^ A m c J^ A m + C P , 



(5) 



where K AC and k AC are the dissociation constant and the 
catalytic constant of active CheA, respectively. 

The dephosphorylation of phosphorylated CheY by 
CheZ (represented by F c in Figure 1C and by Z in the 
equation below) is given by, 
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(6) 



where K Z c and k Z c are the dissociation constant and the 
catalytic constant of CheZ, respectively. 

The dynamics of these processes can be described by 
using a set of differential equations that model the 
NFBLB network depicted in Figure 1C [33]: 
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(9) 



where Oi and 0 2 represent two steady-state values, re- 
spectively corresponding to the two input values I 2 and 
I 2 (Ii = 0.5 and I 2 =0.6 following [33]), and O p is the peak 
value of a transient pulse in response to the input 
change (see Figure 1A). 

Sampling parameter values and numerical simulations 

As in [33], Latin hypercube sampling [35] was used to 
sample uniformly at random the values of kinetic 
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parameters on a logarithmic scale, with the catalytic rate 
constant k being in the range of [10~\ 10 1 ] and the 
Michaelis-Menten constant K in the range of [10~ 3 , 10 2 ]. 
These two parameter ranges were chosen by Ma et al 
[33] to, presumably, minimize computational cost while 
covering previously reported values used to model the E 
coli chemotaxis system [36-38]. In order for our results 
to be comparable with those from previous studies, we 
opted to use the same parameter ranges in this paper. 
For each parameter set, the model (equation (7) for the 
NFBLB model and equation (SI) in the supplemental 
data for the IFFLP model) was numerically simulated 
and those producing the desired adaptation dynamics 
were identified. The original 10 4 parameter sets sampled 
by Ma et al for the NFBLB model produced only eight 
kinetic solutions [33], which are insufficient for discover- 
ing parameter motifs with any statistical significance. To 
remedy this, a 10-fold greater sampling size was used in 
this study; for the IFFLP model, the original 10 4 sam- 
pling produced 131 solutions, enough for the subsequent 
enrichment tests to be carried out. Results from add- 
itional sets of simulation and also from a run increasing 
the sample size to 10 fold, which increased the number 
of kinetic solutions to roughly 10 fold, indicated that the 
parameter motifs reported here had been reliably de- 
duced (see Discussion). 

Following Ma et al. [33], we discarded those parameter 
sets that render the model to produce extremely small 
steady-state values, persistent oscillations, weakly damped 
oscillations, and exceedingly long transient dynamics. For 
each of the remaining parameter sets (46,715 for NFBLB 
and 6,073 for IFFLP), the sensitivity and precision scores 
(i.e. equation (8) and equation (9) respectively) were calcu- 
lated. We said a particular parameter set was a kinetic so- 
lution to the model if its sensitivity score was greater than 
1 and its precision score greater than 10, as these criteria 
have been used to define perfect adaptation [33]. It can be 
shown from equation (8) and (9) that these thresholds 
were chosen to ensure a stimulus of at least 20% of the 
initial steady-state value, and the system can return back 
to this value within an error of 2%, in consistence with 
those used in experimental measurements (Khan et al. 
[39] and Alon et al. [40]). More stringent thresholds re- 
duced the value range of the parameter motifs, but the re- 
sultant minor changes did not significantly affect the main 
findings (see Discussion). We used computer programs from 
Ma et al. [33] and Matlab software (version 7.6.0.324, release 
R2008a) [41] available at http://www.mathworks.com, to im- 
plement a computational pipeline to simulate both the 
NFBLB model and the IFFLP model numerically. We vali- 
dated our simulation pipeline by reproducing the 
numerical results of [33] and also the steady-state solu- 
tion of a four-node transcriptional regulation cascade 
of [42]. 



Enrichment test 

In order to see if there exists an underlying pattern in 
the values of kinetic parameters for the NFBLB model to 
yield perfect adaptation dynamics, we obtained the kin- 
etic solutions and plotted the distribution of parameter 
values for each kinetic parameter. To test whether the 
resulting distributions of parameter values were enriched 
in the sets of kinetic solutions among those 10 5 (NFBLB) 
parameter sets sampled, we adopted an enrichment test 
that has found many uses in genomics sciences [43,44]. 
For a given kinetic parameter, let N be the total number 
of parameter values sampled (here, N = 10 5 for NFBLB 
model), and let y t be the number of those in N belonging 
to the i th value class in the logarithm scale (i= 1, 2, 
5); furthermore, let M be the number of kinetic solutions 
(here, M = 74 for NFBLB model), and x t be the number 
of those in M belonging to the same z th value class. Par- 
ameter values belonging to the 1 st value class are within 
the interval [10~ 3 , 10" 2 ], those belonging to the 2 nd value 
class are within [10~ 2 , 10" 1 ], and so on and so forth, with 
the 5 th value class containing parameter values within 
[10 1 , 10 2 ]. The five value classes of kinetic parameters 
may correspond to varying strengths of enzymatic reac- 
tions that can be measured and classified experimentally. 
Doubling the number of value classes resulted in a simi- 
lar, albeit finer, map of parameter motifs, but did not 
alter the conclusions reached (see Discussion). 

Under the condition of M being sampled independ- 
ently and uniformly at random without replacement 
from N, the probability of observing x t by chance follows 
a hypergeometric distribution [43,44]: 



p(n=x i \y i ,M,N) 



where i = 1, 2, 3, 4 and 5, and 
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As in the enrichment test employed in genomics sci- 
ences [43,44], we can then compute the ^-value to meas- 
ure the statistical significance of the likelihood for 
observing x t when the null distribution (equation (10)) is 
assumed to be the true count distribution. In this study, 
we used a p- value threshold of 10~ 3 to decide the statis- 
tical significance of the enrichment test. 

Thus, for each Michaelis-Menten constant K, we car- 
ried out five independent enrichment tests, each for each 
value class, and for each catalytic rate constant /<, we 
carried out two independent enrichment tests, one for 
the 3 rd value class (i.e. [10"\ 10°]) and the other for the 
4 th value class (i.e. [10°, 10 1 ]), due to the smaller range 
of parameter values for h (see above). If an enrichment 
test was statistically significant (p- value < 10" 3 ), a motif 
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in the form of value class was assigned to the kinetic 
parameter tested. To sum up, the outcome of this ana- 
lysis produced motifs of kinetic parameters, which tells 
us whether a particular kinetic parameter is biased to- 
wards any specific value class on the logarithm scale, or 
none at all 

Functional association test and parameter 
inter-dependence 

For each kinetic parameter, values from all the kinetic 
parameter sets (46,715 for NFBLB and 6,073 for IFFLP) 
were partitioned into two groups. The first group, called 
the motif group, consists of parameter values belonging 
to the biased value class (es) as identified by the above 
mentioned enrichment test on kinetic solutions; and the 
second group, called the non-motif group, comprises 
those parameter values that are not in the motif group. 

We then tested whether the motif group exhibited 
higher sensitivity or precision scores than the non-motif 
group by comparing the score medians of the two 
groups. Since the scores were not normally distributed, 
we used the Mann- Whitney £/-Test [45] for the analysis. 
The test produced a z-score, and we said a particular 
kinetic parameter is positively associated with the sensi- 
tivity or precision parts of the adaptation dynamics if the 
corresponding z-score is greater than 3.29 (i.e. the upper 
boundary of the critical value of the 99.9% confidence 
intervals). In this study our focus was on finding param- 
eters that can significantly improve the function, al- 
though some of the parameters may exhibit a large 
negative z-score indicative of a negative role in the 
function. A bipartite kinetic-functionality network can 



then be constructed to display the associations be- 
tween the kinetic motifs and the functionalities 
identified. 

Finally, we investigated the cooperation between 
every pair of kinetic parameters. Specifically, we took 
the kinetic parameters from those 74 kinetic solutions 
(131 for IFFLP) and performed Pearson correlation test 
between the values of a pair of kinetic parameters that 
appeared in the kinetic motifs identified. We said two 
parameters are correlated if the p-value of the correl- 
ation test is less than 0.05. 



Results 

Detecting kinetic motifs 

As described in the Methods, we found only 74 sets of 
parameters from 10 5 randomly sampled sets that could 
satisfy the criteria of perfect adaptation dynamics for the 
NFBLB model. Figure 2 shows the distributions of these 
74 sets of parameter values. We can see that while some 
parameters, k FBB especially, were limited to be within a 
relatively small range of values, others, like k FC c an d 
K AC , saw a distribution covering nearly the entire range 
of the values sampled. On the whole, catalytic rate con- 
stants kfBB an d k BA were biased towards the 3 rd and the 
4 th value class, respectively, while all five Michaelis- 
Menten constants were biased towards one or more of 
the first three value classes. Other catalytic rate con- 
stants, namely k AB , k AC and k FC c> did not show apparent 
biases towards any vale classes. These observations were 
quantified with statistical significance by the enrichment 
tests; the results, summarized in Table 1, are consistent 
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Figure 2 Distributions of kinetic parameters. Distributions of parameter values obtained from a total of 74 kinetic solutions exhibiting perfect 
adaptation for the NFBLB model. Each tick on the x axis is a specific catalytic rate constant k or Michaelis-Menten constant K, and the values of 
the parameters are in power of 10, which are divided into five value classes as indicated at the right of the figure. On each data box, the 
contracted center is the median, while the edges of the box are the 25 th and 75 th percentiles of the distribution. 
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Table 1 Results of parameter enrichment test for the NFBLB model 



Parameter 




1 st class 




2 nd class 
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4 th class 




5 th class 




[10' 3 , 10 


- 2] a 
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k B A 
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1 .0E + 00 
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73 


50,000 


<2.2E-16 
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50,000 


1 .0E + 00 
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50,000 
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Kab 


44 


19,998 


9.4E-1 1 


25 


19,999 


1 .7E-03 


4 


20,000 


1 .0E + 00 


1 


20,000 


1 .0E + 00 
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20,000 


1 .0E + 00 


Kac 


26 


20,000 


7.4E-04 


17 


19,999 


2.1E-01 


17 


20,000 


2.1E-01 


6 


20,000 


1 .0E + 00 


8 


20,000 


9.7E-01 


Kba 


33 


19,994 


4.3E-07 


28 


20,000 


1.1E-04 


11 


20,000 


8.3E-01 
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20,000 


1 .0E + 00 
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1 .0E + 00 
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17 
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2.1E-01 
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1 .0E + 00 
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20,000 


1 .0E + 00 


0 


20,000 


1 .0E + 00 
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a ln square brackets are the intervals of parameter values for the indicated class. 

b Out of M (=74) kinetic solutions, x is the number of solutions with the value of the indicated parameter belonging to the indicated value class. 
c Out of a total of N (=10 5 ) parameter sets sampled, y is the number of sets with the value of the indicated parameter belonging to the indicated value class. 
d An enrichment test is considered statistically significant if its p-value < 10~ 3 and is highlighted in boldface. Based on p-values (using 10~ 3 as threshold), an 
enrichment state, i.e. motif, was assigned. 



with the distributions of parameter values shown in 
Figure 2. 

Functional roles of kinetic motifs 

For each of the seven kinetic parameters showing bias to 
at least one value class with statistical significance as de- 
termined in Table 1, we carried out Mann- Whitney 
U-test [45] to find out whether a motif (i.e. a preferred 
value class) of that parameter is involved in the sensitivity 
function of the adaptation dynamics, or, correspondingly, 
whether the median of the sensitivity scores for the motif 
group is significantly larger than that for the non-motif 
group (see Methods). The same analysis was then repeated 
for precision. 

The results, summarized in Table 2, indicate that three 
kinetic parameters, K AC , Kfbb and I<fcc> were highly 



significant in improving sensitivity scores, and six kinetic 
parameters, k BA , k FBB , Kab, K ac , K ba and K FC c> were highly 
significant in improving precision scores. Furthermore, for 
the majority of parameters, their motifs were either associ- 
ated with improving sensitivity or precision scores, but not 
both. Two interesting exceptions are K AC and K F cc as their 
motifs tended to reflect improvements in both functional- 
ities (Table 2). These observations are generally in accord 
with the results from an analytical analysis of the rate 
equations (see below and supplemental data). Note that 
sensitivity and precision are two conflicting dynamic pro- 
cesses with the former requiring the system to deviate 
from steady state abruptly and the latter requiring the sys- 
tem to return to the original steady state in a timely fash- 
ion (Figure 1A). Therefore, theoretically speaking, a 
parameter that improves one function will also have a 



Table 2 Functional association test results for the NFBLB model 



Parameter 


Motif [m] 


Non-motif [~m] 




Precision test b 






Sensitivity test b 
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z d 
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z d 
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[10°,10 1 ] 


[10" 1 ,10°] 


1.81 


1.35 


30.88 


-1.18 


-0.76 


-8.75 


PR 


kFBB 


[10" 1 ,10°] 


[10°,10 1 ] 


1.68 


1.45 


16.71 


-1.16 


-0.80 


-6.67 


PR 


Kab 


DO" 3 , 10~ 2 ] 


[10~ 2 10 2 ] 


1.59 


1.42 


6.60 


-1.40 


-0.90 


-6.65 


PR 


Kac 


[10~ 3 , 10~ 2 ] 


[10" 2 10 2 ] 


1.57 


1.52 


9.23 


1.23 


-1.38 


26.57 


PR,SN 


Kba 


DO" 3 , 10" 1 ] 


[10" 1 ,10 2 ] 


1.76 


1.29 


35.34 


-1.41 


-0.68 


-12.72 


PR 


Kfbb 


DO" 3 , 10" 2 ] 


[10" 2 10 2 ] 


1.47 


2.05 


-34.32 


-0.61 


-1.06 


3.95 


SN 


Kfcc 


[10" 1 , 10°] 


DO" 3 , 10" 1 ] and [10° 10 2 ] 


1.63 


1.34 


16.98 


3.30 


-2.31 


28.37 


PR,SN 



a ln square brackets are the intervals of the parameter values, 'm' is the motif group and '~m' the non-motif group. Note that for the kinetic parameters {k AB , k AC , 
and k FC c) showing no apparent bias towards any value classes, the statistical tests were not conducted because their parameter values could not be partitioned 
into motif group and non-motif group (see Methods). 

b "Pr m " ("Sn m ") is the mean logarithm value of precision (sensitivity) scores for the motif group, and "Pr„ m " ("Sn_ m ") is the same but for the non-motif group. 
C "PR" ("SN") indicates that the corresponding kinetic motif is statistically significant (z-score ^3.29) in improving precision (sensitivity). 
d z-score greater than 3.29 (99.9% confidence interval) is highlighted in boldface. 
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negative effect on the other function, as reflected generally 
by the opposite signs of the z-scores obtained from the 
sensitivity and precision tests (Table 2). Finally, the above 
findings can be succinctly captured by drawing a kinetic- 
functionality network (Figure 3), which highlights both the 
functional roles and constraints (i.e. the enriched value 
classes) of kinetic parameters. 

Cooperation of kinetic parameters 

Next, we asked the question of whether kinetic parame- 
ters work independently or in a cooperative way with 
each other when contributing to the system s adaptation 
dynamics. Figure 4 shows the results of correlation tests 
performed on all pairs of the seven kinetic parameters 
that exhibited value class biases. We identified one sig- 
nificant positive correlation (p-value < 0.05) between pa- 
rameters k FB B and K BA , and four significant negative 
correlations in {k FBB , K AB ), {k FBB , K FC c)> (I<ab> I<fbb) and 
(K B a> K F cc) pairs. Interestingly, with the exception of 
I<fbb> most of these correlated parameters contributed 
significantly to systems precision (Figure 3), suggesting 
that they work in a cooperative manner in the precision 
mechanism of adaptation dynamics. In contrast, kinetic 
parameters contributing to the systems sensitivity (i.e. 
Kaci Kfbb and K FC c> see Figure 3) were not correlated 




Sensitivity 



Precision 



-3-2-1012 

Parameter value 
(power of 10) 



Figure 3 A kinetic functionality network. A bipartite network 
connecting kinetic parameters to functionalities (sensitivity and 
precision) of the adaptation dynamics. On the left are kinetic motifs 
emerged from the enrichment tests (see Methods), where filled 
boxes represent enriched values bounded by the indicated power of 
10 for the indicated parameter. On the right are different 
functionalities (sensitivity and precision) of adaptation dynamics. 
A connection between a kinetic parameter and functionality was 
established if the association between the two was determined to 
be significant in the statistical test (see Methods). 

V ) 



with each other, implying that they function independ- 
ently in the sensitivity mechanism. If the mathematical 
model employed can indeed simulate the mechanics of 
adaptation in real biological systems, a corollary of 
these findings is that precision seems to be a more 
complicated mechanism than sensitivity in the system, 
thus the former requires many parameters to work to- 
gether in order to achieve the desired level of high 
precision. 

Experimental data for kinetic parameters of E coli 
chemotaxis 

The NFBLB model is equivalent to the E. coli chemo- 
taxis model with nodes A, B and C corresponding to the 
CheA-bound receptor complex, CheB and CheY, re- 
spectively. Chemotaxis of E. coli has been studied ex- 
perimentally and parameters relating to individual 
processes involved in the model have been estimated 
[36-38,46]. Intriguingly, as described below, the kinetic 
motifs observed here for the NFBLB model compared well 
with the experimental data for the E. coli chemotaxis 
(Table 3). 

To facilitate and simplify the comparison, we dictated 
that parameters with a value in the first three value clas- 
ses (i.e., < 10°) are small-value parameters, while those 
in the 4 th and 5 th value class (i.e. ^10°) are large-value 
parameters. Note that in enzyme kinetics, an enzyme 
with a large K value (Michaelis-Menten constant) indi- 
cates a weak binding affinity to its substrate, while a 
large k value (catalytic rate constant) implies the occur- 
rence of a rapid catalytic event [49]. For sensitivity, the 
empirical estimates for K AC and Kfcc> a fter being nor- 
malized with the concentration of CheY, were 0.36 and 
0.006 respectively (Table 3). This echoes our finding that 
both parameters should take small values. K AC and K FC c 
are respectively involved in the rates of phosphorylation 
(i.e. activation) and dephosphorylation (i.e. deactivation) 
of the response regulator CheY. Goldbeter and Koshland 
[50] explored a simple model of enzyme reaction and 
found that if the activating and deactivating enzymes op- 
erate at saturation where the substrate concentration 
does not affect reaction rate, then an ultra- sensitivity re- 
sponse is observed. Similar arguments may apply to 
chemotaxis in E. coli: namely, in order to produce sensi- 
tivity dynamics, both phosphorylating and dephosphory- 
lating agents (i.e. CheA and CheZ, respectively) must be 
saturated, implying a situation where the concentration 
of the substrate (i.e. CheY) cannot alter the reaction rate. 
A close inspection of the rate reaction for node C in 
equation (7) suggests that both K AC and K FC c taking on 
small values can fulfil this condition. Note that if both 
K AC an d K F cc are small, then C> >K FC c and (1-C)> > 
K AC such that 
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Figure 4 Correlation between kinetic parameters. The Pearson correlation coefficients between pairs of the seven parameters that exhibited 
value class enrichment are shown in the top-right triangle, where a box is colored black if the corresponding correlation is significant 
(p-value < 0.05). In the bottom left triangle, the scatter plots of the paired parameters are shown. On the diagonal are occurrence distributions 
of individual kinetic parameters. 



Table 3 Experimental data for the kinetic parameters of E. coli chemotaxis 



Description 






Reaction 3 
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Notes c [Reference] 
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1.2 


0.08 


k exp = 1 .2 s" 1 [46] and if xp = 0.39 juM [46]; k BA = 
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0.281 


k ex P = 3 2 s -i [37] and ^xp = 1 _ 405 ^ M [37]; ^ 
and K AB = K exp /[B t ]. 


= k ex P 


CheB 


B p 




► B 




V B 2 


0.35 




k exp = 0.7 s" 1 at 35°C, or 0.35 s" 1 at 25°C [47]; k F 
k exp . 


~BB ~ 


Dephosphotransfer 


















CheY 

Phosphotransfer 


A m 
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K AC ^ \ m (~ ^ AC > /\ 




V C1 


650 


0.36 


k exp = 650 s" 1 [48] and if xp = 6.5 fdM [48]; k AC = 
and K AC = k exp /[C t l 


\e xp 


CheY 

Dephosphotransfer 


Fc 


+ c p < 


^ F C C P ^^ 


Fc + C 


Vc2 


30 


0.006 


k exp = 650 s" 1 [34] and /f xp = 0.1 jjhA [34]; k FCC = 
and K FCC = ^ xp /[C t l 


k exp 



a Superscript m denotes methylated form and superscript p denotes phosphorylated form. 

b Each biochemical reaction is equivalent to the process of activation (v n7 ) or deactivation {v n2 ) {n-A, B, and O as indicated in the NFBLB model (see equation (7) 
in Methods). 

c Total concentrations for the CheA-bound receptor complex ([A t ]), CheB ([B t ]) and CheY ([C t ]) are 5.0 uM [34], 2.27 uM [37] and 17.9 uM [34], respectively. 
Superscript exp denotes experimental measurement. 
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and the concentration of node C (i.e. CheY) disappears 
from the rate equation completely. 

The kinetic parameter K FBB is involved in the deactiva- 
tion of node B in the NFBLB model and K FBB is also 
found to be biased to small values (1 st value class motif 
in Table 1) for the model to exhibit sensitivity dynamics. 
Although node B of the NBFLP model corresponds to 
CheB in the chemotaxis of E. coli, to the best of our 
knowledge there are no natural enzymes reported to de- 
activate CheB p in the manner suggested by the NFBLB 
model. Perhaps an unrevealed enzyme exists to dephos- 
phorylate CheB p , the phosphorylate form of CheB, or 
there might be other mechanism of CheB p dephosphory- 
lation that is more complicated than the current know- 
ledge can offer. 

As for precision, the empirical estimates for K AB (nor- 
malized with the concentration of CheB), k FBB , K BA (nor- 
malized with the concentration of receptor complex) 
and k BA are 0.281, 0.35, 0.08 and 1.2, respectively 
(Table 3). These are also in agreement with our findings 
(Table 1) that K AB , k FBB and K BA are biased towards 
small values, and k BA is constrained to large values. K AB 
and k FBB are involved in the phosphorylation and de- 
phosphorylation of CheB, respectively. According to Ma 
et al. [33], K AB is constrained to small values by the 
topological features of the NFBLB model such that the 
rate equation of node B (i.e. CheB) is independent of 
the input level /; this then implies the system is in a 
stable state independent of the initial perturbation and 
thus is able to maintain high precision level. From equa- 
tion (7), a small value for k FBB can reduce the dephos- 
phorylation rate of CheB, this in turn increases the 
deactivation rate of the receptor complex, thereby ensu- 
ing high precision. Finally, K BA and k BA play a part in 
the demethylation of the receptor complex, and intui- 
tively the rate of its demethylation must be great (e.g. 
with a large k BA and a small K BA ) such that the perturb- 
ation initiated by the input signal can be mitigated in 
order to maintain high precision. 

These observations of key parameters of the NFBLB 
model are furthermore supported by a number of ex- 
perimental findings on the chemotaxis system of E. coli: 
1) Kfcc having a large effect on the sensitivity part of 
the system dynamics is in agreement with CheZ playing 
an important role in adjusting the concentration of CheY 
[43]; 2) K AC being an important kinetic parameter in 



affecting sensitivity agrees well with the receptor com- 
plex being a major contributor to sensitivity [44]; 3) 
CheB mutants being far less sensitive than the wild type 
due to the functional abnormality of the receptor com- 
plex in those mutants [45] implies that K FBB is an im- 
portant parameter affecting sensitivity; 4) that 
phosphorylated CheB can increase the rate of receptor 
demethylation and thus speeds up adaptation [51] sup- 
ports the finding that k BA , k FBB , K AB and K BA were all im- 
portant parameters in contributing to the systems 
precision mechanism. 

Discussion 

We have developed here a methodology using parameter 
profiling and enrichment statistical tests to uncover not 
only the sets of kinetic parameters with which a model 
produces user-specified system dynamics, but also motifs 
(i.e. enriched value classes) of these parameters and their 
associations with specific functional aspects of the sys- 
tems dynamics. For these tasks, conventional methods 
usually focus on identifying the "best" parameter set in 
fitting the empirical data and on using complicated ana- 
lytical explorations of models (e.g., sensitivity analysis) 
[52-54] or by a laborious local approach examining how 
changing the parameters one at a time would affect sys- 
tems dynamics [55-57]. Note that the sensitivity analysis 
of conventional methods (not to confuse with the sensi- 
tivity of adaptation (Equation 8) studied in this work) is 
usually carried out for one specific output dynamics, 
whereas our profiling approach investigated the sensitiv- 
ity of functional elements (sensitivity and precision of 
adaptation) for a collection of outpout dynamics (those 
that qualified as a perfect adaptation; there were 74 for 
the NFBLB model). Interestingly, the kinetic parameters 
showing value class biases (Table 2) are those exhibiting 
non- negligible sensitivity indices obtained from analyt- 
ical derivation (Additional file 1: Equation S2) or from 
the numerical method implemented in AMIGO [58] 
(Additional file 1: Figure SI), the latter two being com- 
puted using the output dynamics and the parameter set 
of a randomly chosen one of the 74 kinetic solutions. 
Note that bootstrap-derived distributions of the kinetic 
parameters and their confidence intervals for this par- 
ticular kinetic solution were generally in accord with the 
kinetic motifs deduced from the 74 kinetic solutions 
(Additional file 1: Figure S2). In summary, by profiling 
the parameters as a whole our method takes a global 
view to find not just one but clusters of viable parameter 
sets, thus moving a step further to account for the com- 
plexity of biological systems. Although a number of 
global-view approaches have recently been developed to 
sample from large and high dimensional parameter 
spaces, including a combined global and local explor- 
ation [15] and an approach with model checking on 
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partitioned regions of parameters [17], these studies do 
not make inferences on whether motifs of parameters 
exist and how they might contribute to specific elements 
of the system dynamics. Our approach therefore offers a 
simple framework to systematically characterize and elu- 
cidate the functional contribution of kinetic parameters 
in a biological network. 

The kinetic motifs obtained are quite robust in that 
only minor differences in their resolution were resulted 
from independent sampling runs (Additional file 1: 
Figure S3), different thresholds of the adaptation objec- 
tives (Additional file 1: Figure S4), and added number of 
value classes (Additional file 1: Figure S5). Further ana- 
lysis showed that, for NFBLB and to satisfy the statistical 
significance of p-value < 10" 3 , 8 x 10 4 (rather than 10 5 ) 
sampling runs were required to converge and stabilize 
the kinetic motifs (Additional file 1: Figure S6). The 
computational cost of our method is dominated by the 
sampling and simulation step (see an analysis Additional 
file 1: Section III). To investigate the difficulties that will 
inevitably arise from larger networks for our method, we 
artificially linked two modules of NFBLB together (see 
Additional file 1: section IV) while requiring the system 
to produce the same adaptation dynamics as before. The 
results showed that 5 times of sampling/simulations 
(4 x 10 5 vs. 8 x 10 4 ) were needed to produce stabilized 
kinetic motifs (Additional file 1: Figure S7) for the twice- 
sized network (20 kinetic parameters vs. 10). However, it 
should be noted that the required number of sampling/ 
simulations depends on many factors, including the spe- 
cified output dynamics, the level of statistical signifi- 
cance desired, and the network topology (e.g., the IFFLP 
model needed 10 times less number of sampling/simulations 
than the NFBLB model to exhibit stabilized kinetic motifs, 
despite they both having 10 kinetic parameters). Note 
that our approach is quite general and can be inte- 
grated with other approaches. For instance, the Latin 
hypercube sampling (LHS) used here (and as in [33]) 
can be replaced by other methods, such as those men- 
tioned above, to identify the kinetic solutions needed in 
the subsequent enrichment tests. Results from an ex- 
periment of combining LHS and genetic algorithm 
(GA) showed that 10 3 of LHS sampling followed by 100 
generations of GA optimization yielded a similar set of 
kinetic motifs for NFBLB (Additional file 1: Table SI), 
suggesting that optimization can help to find more so- 
lutions from a smaller initial parameter set, but in this 
case the hybrid approach did not reduce the computa- 
tional cost (since it also needed a total of 10 5 model 
evaluations), and could miss some of the marginally sig- 
nificant motifs (Additional file 1: Table SI). While fur- 
ther research is required to fully address the 
dimensionality problem of scaling up the system, com- 
plex biological networks are known to be composed of 



simple recurrent structural components [59-61]; a dee- 
per understanding of these network components is an 
important first step toward a better understanding of 
the assemblage and functioning of a much larger 
system. 

One interesting observation from the case studies of 
the NFBLB and the IFFLP (presented in Additional file 
1: section VI) adaptation model is that the majority of 
parameters seem to contribute to only one single func- 
tionality (i.e. either sensitivity or precision, see Figure 3, 
Additional file 1: Figure S10). If the mathematical 
models are indeed mechanistic, this could have import- 
ant biological implications as follows. The behavior or 
dynamics of a biological system is likely to be the result 
of intricate interactions of many biological processes. If 
a biological process has a major influence on all aspects 
of the systems behavior, then any changes to such a 
process may have a drastic impact on the system. Modu- 
larity is ubiquitous in biological systems as it provides an 
effective mechanism to corral damaging perturbations to 
local consequences [62]. Thus, to ensure system robust- 
ness, evolution might have favored a biological system 
with a fine division of labors (i.e. modularity) among dif- 
ferent biological components or processes [63-65]. Here, 
in the kinetic-functionality association, we may have un- 
covered yet another example of nature s modularity design 
manifestated in the organization of kinetic parameters. 

Previous studies have shown that certain types and ar- 
rangements of network structures are required to pro- 
duce certain types of system dynamics [33,66-68]; here, 
we have shown that, for a given network structure, cer- 
tain types of values, or motifs, also exist for kinetic pa- 
rameters in order to achieve specific system dynamics. 
Our results suggest that, as has been noted by others 
[69,70], system dynamics can put constraints on the 
values of kinetic parameters. The discovery of these mo- 
tifs underscores the intricate inter-relationships be- 
tween structure (i.e. topology) of the biological 
network, kinetic parameters of the reactions involved, 
and the function of the biological system. Delineation 
of these relationships by methods such as the one de- 
veloped here, which is general and can be applied to 
other types of well defined dynamics, will greatly ad- 
vance our understanding of the design principles of 
prototype biological systems. 

Conclusions 

An increasing number of studies have revealed that 
complicated biological systems often share simple and 
universal design principles that are more understandable 
to biologists. The identification of motifs in biological 
networks is a prime example relating recurrent network 
structures to biological functions. Many studies have 
also argued for the importance of kinetic parameters in 
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determining the dynamics of biological networks, but 
dissecting the association between system dynamics and 
kinetic parameters has been difficult. In this study, we 
have developed a methodology, akin to the enrichment 
analysis of gene expression profiles, to determine 
whether a preference of kinetic parameters adopting cer- 
tain parameter values exists. Such preferences, or kinetic 
motifs, encapsulate the possible roles and functional 
constraints of kinetic parameters. Our analysis on 
models for the adaptation dynamics of the chemotaxis 
of Escherichia coli showed that design principles also 
exist from the perspective of kinetic parameters. Our 
methodology, owning to its generality and simplicity, 
provides a computational framework for understanding 
the kinetic mechanics of prototype biological networks. 

Additional file 



Additional file 1: (I) Mechanisms of functional associations for the 
NFBLB model; (II) Additional analysis for the NFBLB model. Figure 

S1-S6; (III) An analysis of computational cost; (IV) Simulation for a model 
of two linked NFBLPs (NFBLP 2 ). Figure S7; (V) Results for GA-augmented 
simulations of the NFBLP model. Table SI; (VI) Results for the IFFLP 
model Figure S8-S11 and Table S2-S3. 
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